#!/bin/zsh -f

#  set -x  # uncomment to de-bug


version="0.5"

##  mos2fobs.zsh

date_time_stamp="$( command date +'%Y-%m-%d-%H%M' )"



print "#####################################################################################################"
print ""
print "  This script will process the mtz file output from mosflm using CCP4 programs"
print "  including sort, scala, truncate, uniquify, and (optionally) mtz2various to"
print "  produce reduced data in CCP4 format and (optionally) CNS format.  The script employs"
print "  OS X interface dialogs and is meant to be run interactively.  Use this by issuing"
print "  mos2fobs.zsh and optionally the mtz file to be processed.  It requires CCP4 to be installed"
print "  (preferably via fink) and the Camino web browser (the Cocoa version of Mozilla/Netscape)."
print "  It will also make use of phenix.xtriage if it is present."
print ""
print "#####################################################################################################"
print ""

# open -a X11

# Make sure CCP4 is initialized

if [[ -z $CCP4_BIN ]]; then
	if [[ -f $SWPREFIX/bin/init.sh ]]; then
      . $SWPREFIX/bin/init.sh
	fi
fi

if [[ -z $CCP4_BIN ]]; then
	print "I can't find CCP4."
	print "Terminating program.  Make sure ccp4 is initialized."
	return 1
fi

 

# Script to take the mtz output from mosflm and sort, scale, and truncate the data.


# Now define some functions to get needed input parameters from the user

################################################################################
# function GetProjectName lets user specify the CCP4 Project Name (optional)

function GetProjectName {
    osascript << eof-1
        tell app "Finder"
            activate
            set prompt to "Enter a CCP4 PROJECTNAME without spaces (eg: Airborne_HIV) or hit OK to leave blank: "
            set dialogResult to display dialog prompt default answer ""
        end tell

eof-1

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}

################################################################################
# function GetNumAA lets user specify an estimated number of amino acids and GetNumNT
# the number of nucleotides in the
# asymmetric unit of the crystal for quasi-absolute scaling with truncate
# The number isn't critical; I let it default to 150.

function GetNumAA {
    osascript << eof-7
        tell app "Finder"
            activate
            set prompt to "Number of amino acids per MONOMER (eg: 150 for protein; 0 for RNA): "
            set dialogResult to display dialog prompt default answer "75"
        end tell

eof-7

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}

function GetNumNT {
    osascript << eof-71
        tell app "Finder"
            activate
            set prompt to "Number of nucleotides per MONOMER (eg: 25 RNA/DNA bases): "
            set dialogResult to display dialog prompt default answer "25"
        end tell

eof-71

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}

function GetNumASU {
    osascript << eof-72
        tell app "Finder"
            activate
            set prompt to "Guess the number of copies of the above per ASU: "
            set dialogResult to display dialog prompt default answer "1"
        end tell

eof-72

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}

################################################################################



################################################################################
# function GiveLabel lets user customize the labels for F, Sigma, etc.

function GetLabel {
    osascript << eof-5
        tell app "Finder"
            activate
            set prompt to "Enter a subscript name for your dataset (eg: foo  for F_foo): "
            set dialogResult to display dialog prompt default answer "foo"
        end tell

eof-5

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}

################################################################################
# function ChooseFile allows picking from filtered list of files in $PWD
# returns name of chosen file as a string

function ChooseMtzFile {
#
# Change this first line for particular filtering needs:
#
# =========>
#
	filelist=($( command ls -1d *.mtz )) 
#
	if [[ -z $filelist ]];then
		return 1
	fi

	item_list="" 

	for item in "${filelist[@]}" 
	do 
		item_list="$item_list""\"${item}\"," 
	done 

	function filepicker {
		osascript << eof 
			tell app "Finder" 
				activate 
				choose from list {${item_list%,}} with prompt "Pick one: " 
			end tell
eof

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

	} 

	SelectedFile=$(filepicker) 

	print "$SelectedFile"

}
################################################################################

# Function for buttons
# User must supply $buttonDialog

function buttonfuncKey { 
osascript <<-eof
	tell app "Finder"
	activate
	set r to display dialog "$buttonDialog" \
       buttons {"Display", "Dismiss"} default button "Display" giving up after 10
    set favoriteApp to button returned of r
    end
eof

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}


function buttonfuncKey2 { 
osascript <<-eof
	tell app "Finder"
	activate
	set r to display dialog "$buttonDialog2" \
       buttons {"YES", "NO"} default button "YES" giving up after 10
    set favoriteApp to button returned of r
    end
eof

    if [[ $? != 0 ]];then
        print "Cancelled."
        return 1
    fi

}

#################################################################################
#

if [[ $# == 1 ]]; then

   mosflm_mtzfile=$1

elif [[ $# == 0 ]];then

   mosflm_mtzfile=$(ChooseMtzFile) 2>/dev/null
   
    if [[ $? != 0 ]];then
        print "WARNING: mtz file not found"
        print ""
        print "You must run this in a directory containing a mosflm mtz file."
        return 1
    fi
   
   	if [[ $mosflm_mtzfile == false ]]; then
	   print "Selection has been cancelled."
	   return 1
    fi

else 

  mosflm_mtzfile=$1
  print "Reducing $mosflm_mtzfile only.  Ignoring other input."
  print ""
fi

#############################################################################

# Get the Project Name

# PROJECTNAME=$(GetProjectName | cut -d: -f2 | cut -d, -f1)
 
    PROJECTNAME=$(GetProjectName | cut -d: -f3)



   	if [[ $PROJECTNAME == "Cancelled." ]]; then
  	   print "Selection has been cancelled."
  	   return 1
    fi 

print "You have chosen $PROJECTNAME as the project name"
 
#############################################################################

# Get the estimated asymmetric unit content.

# nresidue=$(GetNumAA | cut -d: -f2 | cut -d, -f1)
 
    nresidue=$(GetNumAA | cut -d: -f3)
 


   	if [[ $nresidue == "Cancelled." ]]; then
	   print "Selection has been cancelled."
	   return 1
    fi 


 
# Get the estimated asymmetric unit content.

# nresidueNT=$(GetNumNT | cut -d: -f2 | cut -d, -f1)
 
    nresidueNT=$(GetNumNT | cut -d: -f3)
 


   	if [[ $nresidueNT == "Cancelled." ]]; then
	   print "Selection has been cancelled."
	   return 1
    fi 

# Get the estimated asymmetric unit content.

# nASU=$(GetNumASU | cut -d: -f2 | cut -d, -f1)

    nASU=$(GetNumASU | cut -d: -f3)
 


   	if [[ $nASU == "Cancelled." ]]; then
	   print "Selection has been cancelled."
	   return 1
    fi 



sumresid=$(( $nresidue + 3*$nresidueNT ))

sumasu=$(( ${nASU}*${sumresid} ))

print "Estimate asymmetric unit content as $nASU times the sum of $nresidue amino acids and $nresidueNT nucleotides,"
print "which is very roughly equivalent to $sumasu amino acids for purposes of putting the data on"
print "an absolute scale."

#############################################################################

# Label assignment:


# Assign the labels for the data used in truncate

# fobslabel=$(GetLabel | cut -d: -f2 | cut -d, -f1)
 
    fobslabel=$(GetLabel | cut -d: -f3)
 

   	if [[ $fobslabel == "Cancelled." ]]; then
	   print "Selection has been cancelled."
	   return 1
    fi 
 

print "You have chosen the data label subscript $fobslabel"
 
mtz_prefix=$mosflm_mtzfile:r  

if [[ mtz == $mosflm_mtzfile:e ]]; then

	print "You chose file $mosflm_mtzfile"
    print "Using prefix $mtz_prefix"
	print "Running CCP4 Scala to scale the data .................."

else
	print "There are either no mtz files present in this directory: $PWD"
    print "or you have not chosen any of them because you cancelled the selection."
	print "Program terminating..."
	return 1
fi

################################################################################

#  Sorting data:

################################################################################


sortmtz hklin $mosflm_mtzfile hklout /tmp/junk_${mtz_prefix}_sort.mtz  << end-sort > /dev/null
H K L M/ISYM BATCH
end-sort

if [[ ! -e /tmp/junk_${mtz_prefix}_sort.mtz ]]; then
	print "Sorting error.  Program terminating..."
    print "Please check for possible input mtz file corruption."
    print "(Use for example:   mtzdmp $mosflm_mtzfile )"
	return 1
fi


################################################################################

#  Scaling data:

################################################################################

scala hklin  /tmp/junk_${mtz_prefix}_sort.mtz  \
      hklout /tmp/junk_${mtz_prefix}_scala.mtz \
      << end-scala  >| scala_${mtz_prefix}_${date_time_stamp}_logfile.html
run 1 all
intensities partials # we have few fulls 
cycles 40
# resolution 19.75 2.9
anomalous on           # this is a native set 
#sdcorrection NOADJUST full 5.0 0.00 part 7.0 0.00
#sdcorrection 1.0 0.00   # from a previous run
# try it with and without the tails correction: this is with
#scales   rotation spacing 10  bfactor ANISOTROPIC    
# tails
#reject larger            
#reject 2  # reject outliers more than 3sd from mean
exclude eprob 1e-8    # reject very large observations, if probability
                           #    .lt. 10**-8   

end-scala

buttonDialog="Shall I display the Scala log file in the Camino web browser?"

if [[ $(buttonfuncKey) != Dismiss ]];then
   open -a Camino scala_${mtz_prefix}_${date_time_stamp}_logfile.html
   sleep 5
   buttonDialog=""
fi

if [[ ! -e /tmp/junk_${mtz_prefix}_scala.mtz ]]; then
	print "Scaling error.  Program terminating..."
    print "Please check log file and for possible input mtz file corruption."
    print "(Use for example:   mtzdmp /tmp/junk_${mtz_prefix}_sort.mtz  )"
	return 1
fi

################################################################################

#  Truncate data:

################################################################################

truncate hklin /tmp/junk_${mtz_prefix}_scala.mtz \
         hklout /tmp/junk_${mtz_prefix}_trunc.mtz \
         <<EOF-trunc >| truncate_${mtz_prefix}_${date_time_stamp}_log.html                    
title $mtz_prefix dataset
truncate yes
nresidue $sumasu
PNAME $PROJECTNAME
DNAME $mtz_prefix_$fobslabel
labout  F=F_$fobslabel  SIGF=SIGF_$fobslabel DANO=DANO_$fobslabel SIGDANO=SIGDANO_$fobslabel  
EOF-trunc

buttonDialog="Shall I display the Truncate log file in the Camino web browser?"

if [[ $(buttonfuncKey) != Dismiss ]];then
   open -a Camino truncate_${mtz_prefix}_${date_time_stamp}_log.html
   sleep 5
   buttonDialog=""
fi
 

if [[ ! -e /tmp/junk_${mtz_prefix}_trunc.mtz ]]; then
	print "Truncation error.  Program terminating..."
    print "Please check log file and for possible input mtz file corruption."
    print "(Use for example:   mtzdmp /tmp/junk_${mtz_prefix}_scala.mtz  )"
	return 1
fi

################################################################################

#  Uniquify data:

################################################################################

# Complete dataset and add free-R column for 10% of data.

print ""
print "##################################################################"
print ""
print "  Uniqueify the data; flagging 10% of data for Free-R ............."
print ""
print "##################################################################"
print ""

uniqueify -p 0.1 /tmp/junk_${mtz_prefix}_trunc.mtz \
                 ${mtz_prefix}_unique.mtz >| /dev/null

# Now look at the resulting dataset:

buttonDialog="Shall I examine the new dataset with mtzdmp?"

if [[ $(buttonfuncKey) != Dismiss ]];then

	mtzdmp ${mtz_prefix}_unique.mtz >| check_${mtz_prefix}_mtzdmp.html
	open -a Camino check_${mtz_prefix}_mtzdmp.html
    sleep 5
   buttonDialog=""
fi


# /bin/rm -f /tmp/junk_${mtz_prefix}*.mtz

print ""
print "##################################################################"
print "##################################################################"
print ""
print "    Saving final dataset ${mtz_prefix}_unique.mtz in $PWD "
print ""
print "##################################################################"
print "##################################################################"
print ""

buttonDialog="Shall I display the dataset with HKLVIEW?"

if [[ $(buttonfuncKey) != Dismiss ]];then

	hklview ${mtz_prefix}_unique.mtz >| /dev/null &

	open -a X11  # brings hklview window to the front
    sleep 5
   buttonDialog=""
fi

# Now run phenix.xtriage on the data if it is present and in the path

alias phenix.xtriage="" ; unalias phenix.xtriage

if [[ -x $(which phenix.xtriage | head -n 1 ) ]];then
	print "Running xtriage with results saved to phenix_xtriage_${mtz_prefix}.log"
	phenix.xtriage  ${mtz_prefix}_unique.mtz   \
	                scaling.input.xray_data.obs_labels="IMEAN"  \
	                n_residues=$nresidue  n_bases=$nresidueNT  \
	                n_copies_per_asu=$nASU   \
	                log=phenix_xtriage_${mtz_prefix}.log
	
	# Force this upon the user as it is for their own goddamn good
	open -a Camino phenix_xtriage_${mtz_prefix}.log
fi



################################################################################

#  Optional output of dataset in CNS format:

################################################################################


buttonDialog2="Shall I create CNS formatted datasets?"

if [[ $(buttonfuncKey2) != NO ]];then


# Non-anomalous data:
mtz2various HKLIN ${mtz_prefix}_unique.mtz \
            HKLOUT iso_${mtz_prefix}.hkl << EOF-mtz2cns  >| /dev/null
RESOLUTION 100000 0.5
OUTPUT XPLOR
LABIN FP=F_${fobslabel}    SIGFP=SIGF_${fobslabel}   FREE=FreeR_flag
# LABIN FP=DANO_${mtz_prefix}    SIGFP=SIGDANO_${mtz_prefix}  FREE=FreeR_flag        
END
EOF-mtz2cns


buttonDialog2="Shall I create a CNS formatted dataset for Anomalous data?"

if [[ $(buttonfuncKey2) != NO ]];then

# Anomalous data:
mtz2various HKLIN ${mtz_prefix}_unique.mtz \
            HKLOUT dano_${mtz_prefix}.hkl << EOF-mtz2cns >| /dev/null
RESOLUTION 100000 0.5
OUTPUT XPLOR
LABIN FP=F_${fobslabel}    SIGFP=SIGF_${fobslabel}   FREE=FreeR_flag DP=DANO_${fobslabel}    SIGDP=SIGDANO_${fobslabel}     
END
EOF-mtz2cns
fi

# Edit the cns hkl files to change FOBS, SIGMA and TEST to something else


perl -pi -e "s|FOBS|F${fobslabel}|g"      iso_${mtz_prefix}.hkl
perl -pi -e "s|SIGMA|SIG${fobslabel}|g"   iso_${mtz_prefix}.hkl
perl -pi -e "s|TEST|TEST${fobslabel}|g"   iso_${mtz_prefix}.hkl

perl -pi -e "s|FOBS|F${fobslabel}|g"     dano_${mtz_prefix}.hkl
perl -pi -e "s|SIGMA|SIG${fobslabel}|g"  dano_${mtz_prefix}.hkl
perl -pi -e "s|TEST|TEST${fobslabel}|g"  dano_${mtz_prefix}.hkl



print ""
print "##################################################################"
print "##################################################################"
print ""
print "    Saving final CNS datasets as iso_${mtz_prefix}.hkl"
print "    and (if applicable) dano_${mtz_prefix}.hkl in $PWD "
print ""
print "##################################################################"
print "##################################################################"
print ""



fi

buttonDialog=""
buttonDialog2=""
################################################################################



/bin/rm -i /tmp/junk_${mtz_prefix}*.mtz






